C
      SUBROUTINE CUMAG (X,OUT1,OUT2,NB,NEL,NSOUT,MAG,IGCPL,IGCPS)
C
C     THIS ROUTINE REMOVES EARTH ROTATION AND SENSOR DELAY DISTORTIONS
C     AND MAGNIFIES IMAGERY SURROUNDING A GROUND CONTROL POINT USING
C     CUBIC INTERPOLATION.
C
      LOGICAL*1 X(NB,NEL),OUT2(NB,NSOUT)
      DIMENSION OUT1(NB,NSOUT,4)
      INTEGER POINT(4)
C
      AMAG = MAG
      DO 10 I=1,4
   10 POINT(I)=I
C
C     COMPUTE INPUT PIXELS REQUIRED.  NOTE:  NSOUT SHOULD BE 1 + AN EVEN
C     INTEGER MULTIPLE OF MAG TO BE CENTERED ON AN INPUT PIXEL
      INPIX = (NSOUT+3*MAG-1) / MAG
C
C     COMPUTE BEGINNING LINE AND BEGINNING ELEMENT IN OFFSET COORDINATES
      NL1 = IGCPL - (INPIX-1)/2
      NS1 = IGCPS - (INPIX-1)/2
      OFFNS1 = NS1 - DELAY(IGCPL)
C
C     SELECT REGIONS CENTERED ABOUT GROUND CONTROL COORDINATES
      NL = NL1 - 1
      NREC = 0
   12 CONTINUE
      NL = NL + 1
      READ (10'NL) X
      NROW = POINT(1)
      DO 20 I=1,4
   20 POINT(I) = MOD(POINT(I),4) + 1
C
C     COMPUTE DELAY FOR THIS LINE AND BEGINNING ELEMENT IN CCT COORDS.
      CCTNS1 = OFFNS1 + DELAY(NL)
C
C     FIND BEGINNING CCT PIXEL NUMBER AND DISTANCE TO FRACTIONAL PIXEL
C     AT STARTING PIXEL AND AT ALL SUCCEEDING PIXELS
      DO 50 IB=1,NB
      NS = CCTNS1
      D1 = CCTNS1 - NS
      D2 = AMOD (D1, 1.0/AMAG)
      NPIX = 0
C
   30 CONTINUE
      I1 = X(IB,NS)
      I2 = X(IB,NS+1)
      I3 = X(IB,NS+2)
      I4 = X(IB,NS+3)
      A0 = I2
      A1 = I3 - I1
      A2 = I3 - I4 + 2*(I1-I2)
      A3 = I4 - I3 + I2 - I1
C
C  CUBIC INTERPOLATE A LINE, TO ENLARGE "MAG" TIMES
      INTS = 0
   40 CONTINUE
      D = D1 + INTS/AMAG
      IF (D.GT.1.0) GO TO 35
C
      NPIX = NPIX + 1
      OUT1(IB,NPIX,NROW) = D* (D* (D*A3 + A2) + A1) + A0
      IF (NPIX.EQ.NSOUT) GO TO 50
      INTS = INTS + 1
      GO TO 40
C
C     ADJUST NS AND D1 FOR NEXT CCT PIXEL
   35 NS = NS + 1
      D1 = D2
      GO TO 30
   50 CONTINUE
C
C  ADD "MAG" INTERPOLATED LINES, AFTER OBTAINING 4 INTERPOLATED LINES
C  IN ARRAY 'OUT1'.
      IF (NL-NL1.LT.3) GO TO 12
      DO 70 INTL=1,MAG
      D = (INTL-1)/AMAG
C
      DO 60 IB=1,NB
      DO 60 NP=1,NSOUT
      R1=OUT1(IB,NP,POINT(1))
      R2=OUT1(IB,NP,POINT(2))
      R3=OUT1(IB,NP,POINT(3))
      R4=OUT1(IB,NP,POINT(4))
      CUPIX = D* (D* (D* (R4-R3+R2-R1) + (R3-R4-2.0*R2 + 2.0*R1)) +
     .(R3-R1)) + R2
      CUPIX = AMAX1 (CUPIX, 0.0)
      CUPIX = AMIN1 (CUPIX, 255.0)
      OUT2(IB,NP)=CUPIX+0.5
   60 CONTINUE
C
      WRITE (11) OUT2
      NREC = NREC + 1
      IF (NREC.EQ.NSOUT) RETURN
   70 CONTINUE
      GO TO 12
      END
C
      FUNCTION DELAY (LINE)
C
C     COMPUTE ROTATIONAL AND SENSOR DELAY RELATIVE TO FIRST SWATH
C
C     DEGCEN - LATITUDE AT THE CENTER OF THE LANDSAT SCENE
C     PIXDLY - EQUATORIAL EARTH ROTATION PER SWATH IN PIXELS
C     DEGPLN - CHANGE IN LATITUDE PER SCAN LINE
C     RADDEG - RADIANS PER DEGREE
C     SDELAY - SENSOR SAMPLING INTERVAL BETWEEN LINES (2) / TOTAL NUMBER
C              OF DETECTORS (25)
C     LINESW - LINE NUMBER IN THE SWATH (1 - 6)
C
      DIMENSION SDELAY(6)
      INTEGER CCTLIN, SWATH
      LOGICAL CCT
      COMMON /LANDST/ CCT, LLC, DEGCEN, SAMPOF, LINOFF, AMPL, PHASE
      DATA PIXDLY, DEGPLN, RADDEG, SDELAY /0.6, 0.00071086, 0.01745329,
     .0.0, 0.08, 0.16, 0.24, 0.32, 0.40/
C
      IF (CCT) GO TO 10
      DELAY=0.0
      RETURN
C
   10 CONTINUE
      CCTLIN = LINE + LINOFF
      DEGLAT = DEGCEN + (1170.5-CCTLIN) * DEGPLN
      RDELAY = PIXDLY * COS (DEGLAT*RADDEG)
      SWATH = (CCTLIN-1)/6 + 1
      LINESW = MOD(CCTLIN-1,6) + 1
      DELAY = RDELAY*(SWATH-1) - SDELAY(LINESW)
      RETURN
      END
